﻿using MathNet.Numerics;
using System;
using System.Collections.Generic;
using System.Linq;
using System.Security.Cryptography.X509Certificates;
using System.Text;
using System.Threading.Tasks;
using static TESalgorithm.Constant;
namespace TESalgorithm
{
    public class TES
    {
        /// <summary>
        /// TES算法
        /// </summary>
        /// <param name="BT">五个波段的亮温</param>
        /// <param name="emmisivity">输出的发射率</param>
        /// <param name="Tguess">初始温度</param>
        /// <param name="downatmT">大气下行亮温</param>
        /// <returns>输出地物的温度</returns>
        /// 
        public static double TESalgorithm(List<double> BT,out List<double> emmisivity,double Tguess,List<double> downatmT)
        {

            #region 初始化和初值设置
            //R_b代表去掉大气辐射的辐射亮度
            double[] R_b = new double[5];
            //R代表未去掉大气辐射的辐射亮度
            double[] R = new double[5];
            //大气下行辐射的辐射亮度
            double[] Latm = new double[5];
            //五个波段的猜测温度
            double[] T = new double[5];
            //初始最大发射率设置为0.99
            emmisivity = new List<double> { MaxEmissivity, MaxEmissivity, MaxEmissivity, MaxEmissivity, MaxEmissivity };
            //MMD模块的温度
            double T_MMD = 0;
            //NEM模块的温度
            double T_NEM=0;
            //如果温度的最大值比Tguess要大就替换它
            double Tmax = BT.Max();
            bool NEMstart = true;
            if (Tguess < Tmax) Tguess = Tmax;
            for(int i=0;i<5;i++)
            {
                //根据每个波段对应发射率，求出黑体温度
                R[i] = Equation.radianceOfBand(BT[i], AsterCenter[i], AsterFWHM[i]);
                //大气辐射
                Latm[i] = Equation.radianceOfBand(downatmT[i], AsterCenter[i], AsterFWHM[i]);
            }
            #endregion
            
            #region NEM模块
            //这里不设置终止条件，迭代12次
            for (int iter = 0; iter < NEM_N; iter++)
            {
                for (int i = 0; i < 5; i++)
                {
                    //减去大气辐射的反射
                    R_b[i] = R[i] - (1 - emmisivity[i]) * Latm[i];
                    //估计物体温度
                    T[i] = Equation.Plunk_inverse(AsterCenter[i], R_b[i] / AsterFWHM[i], emmisivity[i]);
                }
                if (NEMstart)
                {
                    T_NEM = Tguess;
                    NEMstart=false;
                }
                else
                {
                    T_NEM = T.Max();
                }
                
                //迭代求解每个波段的发射率
                for (int i = 0; i < 5; i++)
                {
                    double tempRadiance1 = Equation.radianceOfBand(T_NEM, AsterCenter[i], AsterFWHM[i]);

                    emmisivity[i] = R_b[i] / tempRadiance1;
                    //如果发射率超过0.99，就改为0.99
                    emmisivity[i] = Math.Min(0.99, emmisivity[i]);
                }


                //该部分最后得到T_NEM和五个波段的发射率错

                #endregion


                #region Ratio_MMD模块


                //迭代固定次数
                for (int j = 0; j < 5; j++)
                {
                    //发射率比值
                    var beita = Equation.RAT(emmisivity, out double MMD);
                    //经验公式求最小发射率
                    double minEmmisivity =Math.Min(Equation.min_Emmisivity(MMD),0.983);

                    //if (MMD < 0.032)
                    //{
                    //    minEmmisivity = 0.983;
                    //}

                    //得到最小发射率之后，换算成温度
                    //找到最小发射率对应的波段
                    var emMin = emmisivity.Min();
                    int minindex = emmisivity.IndexOf(emMin);
                    //根据最小发射率的比值来迭代求得MMD发射率
                    emmisivity = emmisivity.Select((x,i) => Math.Min(0.99, beita[i]/beita[minindex] * minEmmisivity)).ToList();
                    //分别求五个波段的温度，这里使用的是和NEM一样的方法
                    for (int band = 0; band < 5; band++)
                    {
                        //减去大气辐射的反射
                        R_b[band] = R[band] - (1 - emmisivity[band]) * Latm[band];
                        //估计物体温度
                        T[band] = Equation.Plunk_inverse(AsterCenter[band], R_b[band] / AsterFWHM[band],
                            emmisivity[band]);
                    }
                    T_MMD = T.Max();

                    //迭代求解每个波段的发射率
                    //for (int band = 0; band < 5; band++)
                    //{
                    //    double tempRadiance = Equation.radianceOfBand(T_MMD, AsterCenter[band], AsterFWHM[band]);
                    //    emmisivity[band] = R_b[band] / tempRadiance;
                    //    //如果发射率超过0.99，就改为0.99
                    //    emmisivity[band] = Math.Min(0.99, emmisivity[band]);
                    //    double Rmin_MMD = R[minindex] - (1 - minEmmisivity) * Latm[minindex];
                    //    //得到温度
                    //    T_MMD = Equation.Plunk_inverse(AsterCenter[minindex], Rmin_MMD / AsterFWHM[minindex],
                    //        emmisivity[minindex]);
                    //    //分别求得最后输出的发射率公式(辐射-大气)/(黑体温度对应辐射-大气)，由公式R=B(T)e+(1-e)Atmdown得到
                    //    emmisivity = AsterCenter.Select((x, i) =>
                    //            (R[i] - Latm[i]) / (Equation.radianceOfBand(T_MMD, x, AsterFWHM[i]) - Latm[i]))
                    //        .ToList();
                    //    emmisivity = emmisivity.Select(x => Math.Min(0.99, x)).ToList();
                    //}

                }
            }

            #endregion

            return T_MMD;
            }

        /// <summary>
        /// 用aster的五个波段分别的发射率，求宽波段发射率，默认T=280K
        /// </summary>
        /// <param name="emmisivity"></param>
        /// <returns></returns>
        public static double AsterWideband(List<double> emmisivity,double T=280)
        {
            double totalradiance = 0;
            double totalBBradiance = 0;
            for(int i=0;i<5;i++)
            {
                double centre = AsterCenter[i];
                double FWHM = AsterFWHM[i];
                double em = emmisivity[i];
                double BBradiance = Integrate.OnClosedInterval(lambda =>Equation.Plunk(lambda, T),
                centre - FWHM / 2,
                centre + FWHM / 2);
                totalradiance += BBradiance*em;
                totalBBradiance += BBradiance;
            }
            return totalradiance / totalBBradiance;
        }
    }
}
